function [optPort,FVAL,max_return,CVaR ] = CVaRoptimize( Yij, Q, rf, ...
    alpha, w, v, c)

%                  CVAROPTIMIZE 
%   Function find optimal portfolio and with minimum CVaR and VaR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  Inputs:
% Yij: ---  Senario Prices, without risk free instrument, row = stocks
% Q:   ---  Current Price column vector, without risk free
% rf:  ---  risk free rate over period   (default 0.16%)
% alpha:  - percentile for VaR     (default 95%)
% w:   ---  risk level             (default 10%)
% v:   ---  value constraint for single instrument      (default 20%)
% c:   ---  transaction cost/unit  (default 0)
%                  Returns:
% optPort: optimal portfolio weights including risk free
% FVAL: maximum portfolio value at end of period
% max_return: maximum portfolio return rate
% CVaR: minimum CVaR in percentage of initial portfolio value

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                            
%                  IMPORTANT NOTE
%   Small scale test shows that in order to have feasible results
%   v >= 25%; w >=  1%
%   Large scale: 100 Stocks
%   v can be any value, w >= 8%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% input test
if nargin < 3
   rf = 0.16/100; alpha = 0.95; w = 0.1; v = 0.2; c = 0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function Starts
% Adding risk free instrument
% risk free price is just Nominal Amount, does not matter for portfolio
% return calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rf0 = 100;
[n, J] = size(Yij);  % n stocks, J senarios

% Initialize portfolio, all cash first
x0 = [1; zeros(n,1)];


rf_v = rf0*(1+rf)*ones(1,J);  % Expected cash return is Risk Free rate
Yij = [rf_v; Yij];  % append risk free to Senarios
Ey = mean(Yij, 2);  % getting Ey vector for each instrument
Q = [rf0;Q];        % append risk free to price vector

[n, J] = size(Yij);  % update n

% now we have n+1 instruments, J senarios
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IMPORTANT NOTE: first instrument is Riskfree %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Now we prepare for the 'linprog' program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IMPORTANT NOTE: x vector
% X vector is [ x1:xn, KeXi, Z1:Zj, U1+:Un+, U1-:U1-];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generate 'f' for linprog

f = -[Ey; 0; zeros(J,1); zeros(n,1); zeros(n,1)];
ncol = length(f);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating A matrix for inequality constraint   %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constraint (21)
B1 = w*Q'*x0;  % risk amount of original portfolio
A1 = [zeros(1,n), 1, 1/((1-alpha)*J)*ones(1,J), zeros(1,2*n)];

% Constraint (22)
B2 = -Q'*x0*ones(J,1);
A2 = [-Yij',-ones(J,1),-eye(J), zeros(J,2*n)];

% Constraint (23)
B3 = zeros(n,1);
A3 = [diag(Q')-v*repmat(Q',n,1),zeros(n,1),zeros(n,J),zeros(n,2*n)];

A = [A1; A2; A3];
B = [B1; B2; B3];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating Aeq matrix now   %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test
% Weight should be sum to 1
Beq0 = 1;
Aeq0 = [ones(1,n), 0, zeros(1,J), zeros(1,2*n)];
% end of test

% Constraint (24)
% No transaction cost for risk free instrument (first in Q)
transCon = c*Q';
transCon(1) = 0;
Beq1 = Q'*x0;
Aeq1 = [Q', 0, zeros(1,J), transCon, transCon];

% Constraint (25)
Beq2 = x0;
Aeq2 = [eye(n),zeros(n,1),zeros(n,J),-eye(n), eye(n)];

Aeq = [Aeq0; Aeq1; Aeq2];
Beq = [Beq0; Beq1; Beq2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating lower and upper bound    %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

LB = zeros(ncol,1);
UB = [];


[optX,FVAL,EXITFLAG] = linprog(f,A,B,Aeq,Beq,LB,UB);
% Program will termin
assert(EXITFLAG == 1, sprintf('Exited with flag = %d', EXITFLAG));

% Now get maximum expected return, optimal portfolio, and VaR/CVaR
FVAL = -FVAL;
max_return = FVAL/(Q'*x0)-1;
optPort = optX(1:n);
CVaR = optX(n+1)/(Q'*x0);



end

